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FAST, ADAPTIVE, HIGH ORDER ACCURATE DISCRETIZATION OF THE 
LIPPMANN-SCHWINGER EQUATION IN TWO DIMENSIONS 

SIVARAM AMBIKASARAN, CARLOS BORGES, LISE-MARIE IMBERT-GERARD AND LESLIE 

GREENGARD 

Abstract. We present a fast direct solver for two dimensional scattering problems, where an incident wave 
impinges on a penetrable medium with compact support. We represent the scattered field using a volume potential 
whose kernel is the outgoing Green’s function for the exterior domain. Inserting this representation into the gov¬ 
erning partial differential equation, we obtain an integral equation of the Lippmann-Schwinger type. The principal 
contribution here is the development of an automatically adaptive, high-order accurate discretization based on a quad 
tree data structure which provides rapid access to arbitrary elements of the discretized system matrix. This permits 
the straightforward application of state-of-the-art algorithms for constructing compressed versions of the solution 
operator. These solvers typically require 0(A/’^/^) work, where N denotes the number of degrees of freedom. We 
demonstrate the performance of the method for a variety of problems in both the low and high frequency regimes. 

Key words. Acoustic scattering, electromagnetic scattering, penetrable media, fast direct solver, integral 
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1. Introduction. The problem of acoustic or electromagnetic scattering from penetrable me¬ 
dia arises in a variety of applications, from medical imaging to remote sensing, nondestructive 
testing, sonar, radar and geophysics. In the two-dimensional setting, the governing partial differ¬ 
ential equation is the time-harmonic Helmholtz equation 

Au{x)q{x))u{x) = 0^ x G (1.1) 

where u{x) is the total field and k. is the background wavenumber, and the perturbation (or contrast 
function) q{x) has compact support, say in a domain ft. Note that for q{x) > —1 the solution is 
propagating, while for q{x) < —1 it is evanescent. Using the standard language of scattering theory, 
the total field u{x) can be expressed as the sum of a known incident field u^^^{x) and an unknown 
scattered field u^^^^{x). In order to be well-posed, the latter must satisfy the Sommerfeld radiation 
condition ([Ol. 
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The incoming field is assumed to satisfy the homogeneous equation 


(1.3) 


in some neighborhood D containing Q. From eqs. (O) and ( |1.3[ ), it follows that the unknown 
scattered field satisfies 

= -k^q{x)u''^‘=ix ). (1.4) 

While PDE-based methods can be used to discretize ( |1.4| ) directly (see [22 HU |3S] and the 
references therein), we choose to represent the scattered field using a volume potential 


\x) = V['ip]{x) = ( G^{x,y)i;{y)dy, 
Jn 


(1.5) 



where Jp{y) is an unknown density and G^ix^y) is the outgoing free space Green’s function, given 

by 


G^(x,y) = - y\). 


( 1 . 6 ) 


It is well-known that the operator V in eq. (1.5) is bounded as a map from 1/^(0) to H^{D) and 
compact as an operator acting on [19]. Further, it is straightforward to see that, substituting 

the representation (1.5) into eq. (|1.4[), we obtain the Lippmann-Schwinger integral equation 


'ip{x) + K‘^q{x)V[ip]{x) = f{x) 


(1.7) 


where f{x) = —K‘^q{x)u^^^{x). This is an invertible (resonance-free) Fredholm equation of the 
second kind with a weakly singular kernel. 

Remark 1. The Lippmann-Schwinger equation is often used to denote an alternative formu¬ 
lation HI-- 

u{x) - K^v[uq]{x) = , (1.8) 


derived by convolving the original Helmholtz equation with the governing Greenes function 

Gi^{x,y). The equation (1.1) has two advantages. First, one is often interested in gradients of the 
scattered field. This can he done from (1.5), with full precision by quadrature. Using B 


one 


would need to numerically differentiate the computed solution, with the attendant loss of accuracy. 
The formulation B is also slightly easier to work with, since the contrast function q{x) appears 
as a (left) diagonal multiplier of the volume integral operator, whereas it appears inside the volume 
integral in the classical scheme. 


The Lippmann-Schwinger equation poses several numerical challenges. First, it leads to a large, 
dense linear system of equations for f){x). Second, it may involve a complicated contrast function 
q{x), requiring adaptive mesh refinement for effective resolution. Third, it may be ill-conditioned 
due to multiple scattering once the contrast q{x) and n are large. The literature in this area is 
substantial, and we do not seek to review it here. Some relevant prior work on volume integral- 
based methods, fast solvers, and numerical scattering theory includes 13 011311 El llQl HI m 

Our goal in this paper is to develop a fast, adaptive, high-order accurate, direct solver, which 
shares features with many of the algorithms in the papers cited above. We concentrate, in particular, 
on developing a framework that permits the rapid computation of entries in the governing system 
matrix, once the adaptive data structure has been specified. With this capability in hand, it is 
straightforward to make use of modern hierarchical direct solvers. 

There has been relatively little attention paid to the adaptive discretization issue, and fast direct 
solvers for volume integral equations are typically implemented on either uniform Cartesian or polar 
grids. This is a perfectly sensible approach, especially when first developing fast solvers, whether 
based on separation of variables and FFTs or hierarchical direct methods [llIIllillIiEIllMlig. 

In the next section, we describe a high-order adaptive approach to discretization of the un¬ 
known density in B, followed by a detailed explanation of how one can rapidly compute the 
matrix entries of the fully discretized system. We then discuss a fast solver for the Lippmann- 
Schwinger equation using the HODLR method of HE] and present numerical examples illustrating 
the performance of the scheme. 
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2. Discretization. We assume that we are given a square domain D which contains the 
support of a known contrast function q{x). By inspection of the Lippmann-Schwinger equation ( 0 , 
D clearly contains the support of our right-hand side f{x) = —hz‘^q{x)u^^^{x), and therefore the 
unknown 'ip{x) as well. We subdivide D using an adaptive quad-tree, ensuring that q{x) and f{x) 
are well-resolved, and that the number of points per wavelength is greater than or equal to a user- 
specified parameter in each leaf node, on which we impose a p x p grid (based on a uniform grid for 
p = 4 and based on a tensor product Chebyshev grid for p > 4). The unknowns are taken to be the 
values of 'ip at the p x p grids on every leaf node, and we will seek to enforce the integral equation 
at the same nodes, corresponding to a Nystrom discretization of the integral equation. If we let 
denote the vector of function values at all leaf node grid points and we let / denote the vector 
of right-hand side values at those same leaf node grid points, then the discrete version of eq. (1.7) 
takes the form 


= f. ( 2 . 1 ) 

2.1. The adaptive data structure. The domain D is decomposed hierarchically, as follows. 
Grid level 0 is defined to be the domain D itself. Grid level / +1 is obtained by recursive subdivision 
of each box B at level I into four child boxes (Fig. B is referred to as their parent. We allow 

for adaptivity, so that not all boxes at level I are necessarily divided. We will, however, require 
that the tree satisfy a standard restriction - namely, that two leaf nodes which share a boundary 
point must be no more than one refinement level apart. Refinement is controlled by the following 
two criteria. First, given the wavenumber (Helmholtz parameter) /d, a box is subdivided if it is of 
side length and > 27rM, for a user-specified parameter M. This ensures that there are at 
least M points per wavelength in each linear dimension. Second, we ensure that the contrast q{x) 
and the right hand side f{x) are both well-resolved. For this, suppose we are given a leaf node B 
at level 1. We evaluate the functions q{x) and f{x) at tensor product Ghebyshev nodes on B and 
compute the coefficients of the corresponding pth order Ghebyshev approximation, which we will 
denote by P 23 . We then subdivide the box into four child boxes and compute q{x) and f{x) 
at tensor product Ghebyshev nodes on each one. If the computed values in the chid boxes agree 
with P]s to a user-specified tolerance e, we terminate the refinement at box B. Otherwise, we add 
the children to the data structure at level I + 1 and continue until the approximation criterion is 
satisfied. 
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t 

Bo 
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Fig. 2.1. Box B and its four children 


The procedure described above will not, in general, produce a level-restricted tree. It is straight¬ 
forward, however, to add further refinements until that criterion is satisfied as well (Fig. 2.2). 


2.2. The volume integral. We shall assume that there are M leaf nodes in the tree, denoted 
by {Dm}, so that 


D = U 


M 

m=l 


Dm- 
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(a) Adaptive quad tree before level restriction (b) Adaptive quad tree after level restriction 


Fig. 2.2. An adaptive quad-tree and its level-restricted refinement. 


Thus, the volume integral defining the scattered field (1.5) can be written in the form: 

M 


(x) = v['ip]{x) = ^ I / Hfi{K\x - y\)'ip{y)dy 

m=l 


( 2 . 2 ) 


Note that, since we are using an adaptive tree, the various Dm can be at arbitrary levels in the 
spatial hierarchy. To obtain a Nystrom method, it remains to discuss the computation of a high- 
order accurate discretization of V[ip]. For this, we will require some notation. We let the values of 
the unknown density on leaf node Dm be given by tpmj ~ for j = 1 ,... where Xj is 

one of the grid points on Dm (Fig. |2.2[ ). 



Fig. 2.3. Discretization points for a leaf node Dm and two neighboring boxes at different levels of the hierarchy. 
A uniform grid is used for Jfih order accuracy and a tensor-product Chebyshev mesh for higher order accuracy. 

When clear from context, we write the full unknown vector as for j = 1,...,A/^ where 
N = p^M the total number of unknowns. We write qg = q{xj) and fj = f{xj) for the contrast 
function and and the right-hand side, respectively, at the corresponding grid points. 

2.2.1. Polynomial approximation. In order to obtain p^^ order accuracy, we build a pth. 
order polynomial approximation to the density tp on each leaf node B of length Ls centered at 
(Si, S 2 )- For this, we let x = (^i, ^ 2 ) and let 5j(^i, ^ 2 ) denote a suitable basis for polynomials of two 
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variables of degree p — 1. That is, the bj{^ 1 ,^ 2 ) should span {^ 1^2 • 0 < a + 6 < p with a, 6 G N^}. 

pip + 1 ) 

The number of such basis functions is clearly Np = -^-. In practice, we use the simple 

monomials ^ 1^2 fourth order accuracy and Chebyshev polynomials of the form Ta{^i)T}j{^ 2 ) for 
higher order schemes, scaled to the unit box [— 1 / 2 , 1 / 2 ]^. 

Suppose now that we are given a pth order polynomial defining ' 0 (^i,<^ 2 ) on S: 


Ar„ 


'023(^1:6) = 023(0^2 


2=1 


6 -gl 

Lb 


Lb ) 


(2.3) 


For the tensor product grid points Xi = (Ci, 15 ^ 2 , 2 ) lying in S, we define the interpolation matrix 
Q : with entries Qu by 


Qii = bi 


\ Lb ^ Lb ) 


so that 




'023(6,156,2) = yy QiicBji ) ■ 


1=1 


Note that in eq. (2.3), the basis functions are independent of the level of box B in the quad 
tree, since we normalize by the box dimension Lb- If we let ipB ^ 1^^ denote the function values 
at the tensor product grid points in standard ordering, and cb = ( 0 / 3 ( 1 ), 0 / 3 ( 1 ),..., cb{Np)) G 
then 0/3 Qcb- Let us denote by the pseudoinverse of Q, so that the coefficients cb can be 
computed from 0/3 via 


CB 


= QHi 


The cost of computing by means of the singular value decomposition is negligible. Moreover, 
this is done only once and the pseudoinverse is stored. 

Suppose now that we wish to compute an arbitrary matrix entry Vij corresponding to a pth 


order accurate Nystrom discretization of ^[0] in ( 2 . 2 ). We assume that Xj lies in box B and that Xi 
is an arbitrary point in the adaptive quad tree data structure (including possible B itself). Then, 
for pth order accuracy, the contribution of the jth grid point to the ith target point is precisely 


2=1 


L 


H^P{K\xi - y\)bi 


yi -Bi 2/2 - B2 


Lb 


Lb 


dyidy2 , 


(2.4) 


where y = (^ 1 , ^ 2 )- This follows from the fact that the ith column of Ql is a vector in , consisting 
of the contributions (or projections) of the ith grid point onto each of the Np basis functions. Thus, 


the (i, j) entry of the full matrix in the linear system in eq. ( 2 . 1 ) is given by 


^ij — 6ji T ^ ^(^i)L^ji 


(2.5) 


It is worth emphasizing that the amount of work in computing Vij (or Aij) is independent of 
N (the total number of discretization points). Each entry could, for example, be computed by 
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adaptive quadrature on the fly. We wish, however, to be able to evaluate an arbitrary matrix entry 
in only a few floating point operations. The remainder of this section is devoted to an extremely 
efficient method for this task, based on the geometric relations between Xi and Xj. 

Definition 2.1. The colleagues of a box B are boxes at the same level in the tree hierarehy 
whieh share a boundary point with B. (B is eonsidered to be a eolleague of itself.) Note that eaeh 
box has at most 9 eolleagues (Fig. \2.4^ . 
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Fig. 2.4. Colleagues 

The coarse neighbors of B are leaf nodes that are one level higher than that of B and whieh 
share a boundary point with B. Note that there ean be at most 12 eoarse neighbors (Fig. |^.5| ). 
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(a) 


(b) 



Fig. 2.5. Coarse neighbors 


The fine neighbors of B are leaf nodes that are one level below that of B and share a boundary 
point with B. Note that there ean be at most 12 fine neighbors (Fig. 2.6). 

The separated fine neighbors of B are non-neighboring leaf nodes that are one level below 
that of B and whose parent shares a boundary point with B. There are at most 16 separated fine 
neighbors (Fig. \2.1^ . 

All other leaf nodes are separated from B by at least the length Ls and are eonsidered its far 

field. 


2.3. Far-field interactions. The easiest case to consider is when Xi lies in the far field of box 
B. If we let Xi = (Ci, 15 ^ 2 , 2 ) and denote by Oi the polar angle of Xi with respect to the box center 
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2.6. Fine neighbors 
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Fig. 2.7. Separated fine neighbors 


( 61 , 62 ), the Graf addition theorem [HllO] states that 

00 

Hf\K\xi-y\)= ^1.^2-^2)||)J„(/t||(2/i-^ 1 , 2 / 2 (2.6) 

n= —00 

where Oy denotes the polar angle of the point y E B with respect to the box center. Since the 
dimensions of the leaf node are assumed to satisfy kLjs < 8 , we can truncate eq. ( 2 . 6 ) after about 
1/ + 8 terms with an error of approximately 2 “^. It is easy to verify that 


Vij -jYl - (Si,B2)||)e**' 


(2.7) 


n= — L 


where 


Np 

Mn{j) = '^CnlQl^ 
1 = 1 
7 


( 2 . 8 ) 

































and 


Cnl = [ Jn{K\\{yi - Bi,y 2 - ^ 2 ) 11)6 
Jb 


—iW 


ybi 


Vi -Bi 2/2 - B 2 


Lb 


Lb 


dyidy2 . 


(2.9) 


Note that the integrals in (2.9) are smooth and need to be computed only once per level at negligible 


cost. Moreover, Mn{j) also only needs to be computed once per level. Thus, the cost of computing 
a far field matrix entry to fourteen digits of accuracy is essentially that of evaluating the multipole 
expansion (2.7) with L = 45, which can be done about one million times per second on a single 


2.4. Separated fine neighbors. The multipole expansion for B in the previous section is 
slowly converging for the separated fine neighbors (Fig. 2.7). However, it is straightforward to 
create four child boxes from B and to compute the multipole expansions for each of the four 

'yi-Bi 2 / 2 -^ 2 ' 


children from each basis function bi 


Evaluating the corresponding multipole 


Lb ’ Lb 

expansions at each separated fine neighbor target point yields the desired value Vij. 


2.5. Near field interactions. For points in the near field, we rewrite (2.4) in the form: 

Np 

( 2 - 10 ) 


1=1 


where 




- {yi,y2)\)k 


f yi - Bi y2-B2 \ 

\ Lb ’ Lb ) 


dyidy2. 


( 2 . 11 ) 


It will be convenient to rescale variables so that integrals are always carried out over the unit 


box. Thus, instead of (2.11), we write 




J Ho\i^LBri)hi{yi,y2)dyidy2 


( 2 . 12 ) 


where U = [—0.5, 0.5]^ and = \xi/LB — (^i,^ 2 )|- Note that, in this representation, the integral 
depends only on the relative coordinates of the target point Xi to the unit box U. There are 
only finitely many such possibilities, corresponding to the various target locations in the colleagues, 
coarse neighbors, or fine neighbors. Thus, the number of such interactions is fixed. Xi must be one 
of possible locations in each of 9 possible colleagues, 12 possible fine neighbors and 12 possible 
coarse neighbors. Thus, we seek to compute tables of dimension x 9 x x 12 x p^, and 

p^ X 12 X p^, respectively. For this, we use the fact mm that 


H^o'\z) = Mz)+iYo{z) 


with 


Mz) = 

m=0 


(m!)^ 
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^o(^) = ^ (in (I) + 7 ) Mz) + 


2 f 




1 , lAby 


" -(l+o)^^ + (l+o + 7) " 


TT V(l!)^ 2^ (2!)2 

where 7 is Euler’s constant. This permits us to write 

r rp\ 2 p 

. 2 ) 

p=0 p=0 

where 


2 3^ (3!)2 


H, 


f\.L„r) = (9"’ + (^)‘'7og (9 (2.13) 


2 % 2 % (— 

Cp{z) = ap{z) + —gp(z), dp{z) = —ap{z), ap{z) = , gp{z) = (7 + log(2:) - Hp) ap{KLB). 

TT TT [P-J 

Here, Hq = 0, = 1, H 2 = {I -), = (! + - + -), etc. This permits us to write 




+ 




Pmax p 2p 

^Cp(kLb)/ (0 h{yi,y2)dyidy2 

p=0 

Dmax p 2p 

^ dpiKle) (0 log (0 bi (yi, ^ 2 ) clyirfy 2 1 


(2.14) 


for a suitable choice of Pmax • 

The integrals over U in eq. (2.14) are independent of the wavenumber so we may tabulate 
these values to double precision accuracy using adaptive Gaussian quadrature. For kLb < 8 , it 
is easy to verify that Pmax = bO is sufficient to obtain fourteen digits of accuracy. The quantities 
V-^ are then obtained using 0{pmax) floating point operations. The storage requirements for this 
scheme are: 


p X Np X 9 X (2pj;nax H“ 2) 
p X Np X 12 X (2pjnax H“ 2) 
P^ X NpX 12 X (2pmax + 2) 


for colleagues 
for fine neighbors 
for coarse neighbors. 


Unfortunately, for kLj^ > 4, the computation in ( |2.14 ) can be subject to catastrophic cancella¬ 
tion. (The series oscillates in sign and involves intermediate terms of large magnitude.) This loss of 
precision is easily overcome by subdividing B into four children and generating tables for each child 
box. With this scheme, the storage costs remain negligible and the matrix entries Vij corresponding 
to local interactions are computed using only a few hundred floating point operations. 


3. Fast direct solvers. The reason we have focused in this paper on the rapid computation 
of matrix entries is that a new generation of fast solvers has been developed over the last decade or 
so which permits the inversion of equations such as the Lippmann-Schwinger equation in 
work. All that is required is a suitable ordering of the unknowns and access to matrix entries on 
the fly. For the sake of simplicity, we have chosen to use the HODLR scheme of HE]. This solver 
relies on a hierarchical partitioning of the matrix into a sequence of off-diagonal low-rank blocks 
(from which the acronym is derived). More precisely, the off-diagonal blocks are approximated 
using low-rank factorizations to a user-specified precision e. We refer the reader to the original 
papers for details and to liiaiaiiaiTiiTaiTHiEniiiHiiisiiisiiiaEsiEaiMiiiiiiisi for related 
schemes. All of these fast solvers, of course, require sampling at most 0{N^^‘^) matrix entries. For 
PDE-based analogues, see [211|371|29| and the references therein. 
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Several of the methods above have smaller constants than HODLR for volume integral equations 
with singular kernels, but require a slightly more complicated interface. In fact, in the low frequency 
regime, some of the methods above require only 0{N log N) work (see, for example, [HI [29]). We 
postpone such accelerations to a later date. 

4. Numerical results. In this section, we illustrate the performance of our adaptive solver for 
the Lippmann-Schwinger equation. In each example, we use an incoming plane wave propagating 
in the x-direction of the form and calculate the value of the scattered field on a 

square D which contains the support of the contrast function q{x). 

We will make use of both uniform and adaptive grids. As indicated in sectionin the adaptive 
case, we refine the domain D using a quad-tree, halting refinement of a leaf node when it is 
determined to have resolved the contrast function q{x) and the right-hand side f{x) to a user- 
specified tolerance e. In order to be conservative, we then use a tolerance of e x 10“^ in the 
HODLR solver. The number of terms used in various far field expansions is chosen as in [20|. All 
simulations in this paper were carried out using a single core of a 2.5GHz Xeon processor. 

Example 4.1. Radially symmetric contrast functions 

In our first two examples, we consider a radially symmetric contrast function g(x), from which 
it is straightforward to compute the exact solution using separation of variables, the fast Fourier 
transform, and a high-order accurate ODE solver (see, for example, H). Setting the wavenumber 
K = 40, we will use either a Gaussian with q{x) = or a “flat bump” with q{x) = 

0.5erfc(5(|xp — 1)), where erfc is the complementary error function [T1 


2 ^ 

erfc(x) = —= / e~^ dt. 

Jx 


In Fig. |4. 1(a) I we plot the Gaussian contrast, and in Fig. |4.1(b) we plot the real part of the 
total field after solving the scattering problem. In Fig. |4.2(a) we plot contours of the “flat bump” 
contrast function, in Fig. |4.2(b)| we plot the real part of the total field, in Fig. 4.2(c) we show the 
the adaptive discretization of the domain D, and in Fig. 4.2(d)| we plot the contrast function q{x). 




(a) Contrast function 


(b) Real part of total field 


Fig. 4.1. Gaussian contrast 
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(a) Contrast function 



(c) Grid for adaptive solver 


(b) Real part of total field 



(d) Surface plot of contrast function 


Fig. 4.2. “Flat hump'’ contrast 


We let u denote the total field calculated using our solver, we let u denote the total field 
computed using separation of variables as in m- Table |4.1| shows the performance of the solver on 
a uniform grid. 


It is straightforward to check from Table [4T] that the convergence rate is approximately fourth 
order (for a target either within the support of q{x) or in its exterior). The CPU requirements of 
the solver can also be seen to scale approximately as 


We now solve both the Gaussian and flat bump problems using our adaptive refinement strategy. 
Tables [4^ and [T3| summarize the corresponding numerical results. 

As expected, the adaptive method provides significant advantages: for the same accuracy, the 
problem with a Gaussian contrast function is about ten times faster. 

One of the advantages of using the Lippmann-Schwinger formulation over direct discretization 
of the Helmholtz equation is that Green’s function based methods are not subject to the same 
kind of grid dispersion error. To verify this, in Figs. 4.3(a)| and |4.3(b)[ we plot the error for 
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Table 4.1 

Convergence behavior of Lippmann-Schwinger solver on uniform grid for Gaussian contrast. N denotes the 
number of discretization points, e{x) = \u{x) — u{x)\ is the error at the point x. (0.5, 0) is well inside the support of 
the contrast function and (1,0.5) is a point where q{x) has vanished to machine precision. Time is the time required 
by the fast solver. 


N 

e(0.5,0) 

e(l,0.5) 

Time 

256 

0.8925E + 0 

0.4890.B + 0 

0.14 

1024 

0.1809i; + 0 

+ 0 

1.16 

4096 

1.0938i; - 2 

9.6485E’ - 3 

12.93 

16384 

3.1169E-4 

3.2633.B - 4 

104.1 

65536 

1.4300E - 5 

1.2537.B-5 

844.8 

262144 

8.8874E - 7 

6.4485.B - 7 

6869.7 


Table 4.2 

Performance of adaptive solver for Gaussian contrast, e is the tolerance requested from adaptive discretization 
of the data. The other columns are the same as in Table \4^ 


e 

N 

e(0.5,0) 

e(l,0.5) 

Time 

IE-4 

4096 

1.09E-2 

9.65E’ - 3 

12.3 

IE-5 

4864 

7.89E - 4 

1.20.B - 4 

16.2 

IE-6 

6400 

2.99E - 4 

3.33.B - 4 

30.4 

IE-7 

10240 

4.83E - 5 

8.75E’ - 6 

78.7 

IE-8 

16384 

1.13E-5 

5.35E’ - 6 

190.1 

IE-9 

34816 

2.30E - 6 

•^.mE - 7 

650.6 

IE - 10 

70912 

7.14E-7 

3.16.B-7 

2350.7 

IE - 11 

138688 

4.25E - 8 

1.94.B - 8 

7363.4 


the flat bump contrast function as a function of wavenumber k at the points (0.5,2) and (3,0.5), 
respectively. Each curve represents the error for a fixed value of the parameter kLjs on loaf nodes. 

= 4, 2,1 correspond to approximately 6, 12 and 24 points per wavelength. Note that the error 
does not grow with tz as it would with a fourth order discretization of the PDE. 

Example 4.2. Multiple scattering from well-separated Gaussian bumps 

To illustrate the performance of the scheme on a more complicated example, we assume the 
contrast consists of the sum of 20 Gaussian bumps of the form 

qj{x) = 


where the centers {xj} are randomly located in [—1.5,1-5]^ and a = 0.0013. With wavenumber 
/^ = 60, the domain is approximately 30 wavelengths in each linear dimension. We discretize the 
domain using at least 9 points per wavelength and resolve the contrast with tolerance e. In Eig. 
4.4(a)[ we plot the contrast function, and in Eig. 4.4(b)| we plot the real part of the total field 


calculated using N = 895264 discretization points with about five digits of accuracy. Eig. 4.4(c) 
shows the adaptive grid (at a coarser discretization) and Eig. |4. 4(d) [ presents a surface plot of the 
contrast function. Table [4^ summarizes our numerical results. 

Example 4.3. An example from plasma physics 

Radio frequency waves play an important role in magnetic confinement fusion for both heating 
and diagnostic applications. Radio frequency wave propagation in a fusion reactor can be described 
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Table 4.3 

Performance of adaptive solver for flat bump contrast, e is the tolerance requested from adaptive discretization 
of the data and N is the total number of discretization points. e{x) is the error at the point x. (0.5, 2) is well inside 
the support of the contrast function and (3, 0.5) is a point where q{x) has vanished to machine precision. 


e 

N 

e(0.5,2) 

e(3,0.5) 

Time 

IE-3 

10240 

l.SSE - 1 

9.29F; - 1 

61.8 

IE-A 

40192 

1.07F; - 2 

2.27E - 2 

609.2 

lE-b 

134656 

SATE - 4 

SAAE - 4 

5561.7 

IE-7 

596224 

2.07F; - 5 

1.31F;-5 

75861.5 


Table 4.4 

Convergence behavior of scattering from multiple Gaussian bumps, e is the tolerance requested from adaptive 
discretization of the data and N is the number of discretization points. and ^{u{x)) denote the real and 

imaginary parts of the solution at the point x. 


e 

N 

3?(u(0.5,2)) 

Q{u{0.5,2)) 

3?(m(3,0.5)) 

Q{u{3,0-5)) 

IE-7 

155056 

-1.74310F;- 1 

-2.89411F; - 

1 

-7.73067.B - 1 

-4.49033F; - 1 

if;-8 

272128 

-1.74256F;- 1 

-2.89433F; - 

1 

-7.73039E’ - 1 

-4.49080F; - 1 

if;-9 

441520 

-1.74268F;- 1 

-2.89425F; - 

1 

-7.73045E’ - 1 

-4.49084F; - 1 

if; - 10 

895264 

-1.74266F;- 1 

-2.89425F; - 

1 

-7.73044.B - 1 

-4.49085F; - 1 


using a cold plasma model misiiii], where the plasma is considered a fluid, possibly involving 
several species, subject to electric and magnetic forces. Here we consider only electrons, neglecting 
ions because of the mass ratio. Coupling Newton’s law for the fluid motion with Maxwell’s equations, 
driven by a current proportional to the electrons velocity, leads to various propagation modes. These 
are classified by their orientation with respect to the confining magnetic field Bq. The ordinary 
(O) mode corresponds to an electric field parallel to Bq, propagating in a plane orthogonal to Bq. 
This particular mode is used for reflectometry, a non-invasive method to probe the plasma density. 


The O mode propagation in the cold plasma model reduces to the Helmholtz equation (1.1), 


where the contrast q{x) is proportional to the additive inverse of the electron density, which is 
compactly supported as the plasma is surrounded by vacuum. Reflectometry relies on the fact that 
the incoming wave impinging on the plasma can only propagate if the density is smaller than a 
threshold, called the “cut-ofF’. In equation (0, the cut-off is implicitly defined by q{x) = — 1. 
When the wave reaches the cut-off it cannot propagate further and is reflected back from the plasma. 
In other words, the plasma forms an evanescent medium when the density is higher than the cut-off 
density, i.e. where q{x) < —1. 

In Fig. 4.5(a)[ we plot the contrast for an idealized density profile, defined by 


Q{x,y) 


-{ 


0 , 

-0.05) - g{x,y) cos(0.9y), 


if '0 (t, y) < 0.05, 
if '^(t, y) > 0.05, 


where 


'0(t, y) = 1 — {x — 0.15(1 — — C(1 -f 0.3x)^y^, 

/ {x — Xj)‘^{y — yj)‘^\ 

gix,y) = J2 exp ( - ^j , 

j=l ^ ‘ / 
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(a) Error at the point (0.5,2) as a function of k,. 



Fig. 4.3. Error as a function of n, which determines the size of the domain in wavelengths. 
kLq = 4, 2,1 correspond to using approximately 6, 12 or 24 points per wavelength. 


The values 


with C = 0.4987, and the values of a^, Xj, and yj^ for j = 1,... ,5, are given in Table 4.5 We 


set the wavenumber k, = Sb. To solve this problem, we used N = 1060864 points corresponding to 
approximately 25 points per wavelength. The total field calculated by our solver is shown in Figure 
4.5(b)[ As expected, there is no propagation through the evanescent zone and the incident wave is 


reflected back at the cut-off. Convergence analysis using Richardson extrapolation indicates that 
we obtained more than 6 digits of accuracy. 


5. Conclusion. We have presented a high-order method for the Lippmann-Schwinger equation 
that is of particular value when coupled with modern fast, direct solvers. More precisely, we have 
developed algorithms for rapidly accessing arbitrary elements of the system matrix consistent with 
an adaptive Nystrom discretization. Our approach extends naturally to three dimensions and 
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(a) Contrast function for twenty Gaussian bumps (b) Real part of total field 



Fig. 4.4. Multiple scattering from 20 Gaussian humps 


to other governing Green’s functions, although some of the accelerations have been developed 
specifically for the Helmholtz equation and will require modification in other settings. 

Our numerical examples were computed using a fourth order accurate approach, but the exten¬ 
sion to arbitrary order is straightforward. We are currently working on the coupling of our tools with 
fast solvers that are quad-tree or oct-tree based, rather than k-d tree-based as in our HODLR solver. 
These were discussed very briefly in section This will substantially reduce the constant implicit 
in the notation. We are also working on fast discretization of the Lippmann-Schwinger 

equation in three dimensions for both acoustic and electromagnetic applications. 
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Table 4.5 

Parameters aj, Xj, and yj for the eontrast funetion used in the plasma physies example. 



(a) Contrast function q{x) (b) Real part of total field 


Fig. 4.5. RF wave impinging on plasma. 
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